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Abstract 

We study the prediction of the dark matter power spectrum at two-loop order in the Effective Field Theory 
of Large Scale Structures (EFTofLSS) using high precision numerical simulations. In our universe, short 
distance non-linear fluctuations, not under perturbative control, affect long distance fluctuations through 
an effective stress tensor that needs to be parametrized in terms of counterterms that are functions 
of the long distance fluctuating fields. We find that at two-loop order it is necessary to include three 
counterterms: a linear term in the overdensity, <5, a quadratic term, S 2 , and a higher derivative term, d 2 S. 
After the inclusion of these three terms, the EFTofLSS at two-loop order matches simulation data up 
to k ~ 0.34 LMpc -1 at redshift z = 0, up to k ~ 0.55 /iMpc -1 at z = 1, and up to k ~ 1.1/iMpc -1 
at z = 2. At these wavenumbers, the cosmic variance of the simulation is at least as small as 10 -3 , 
providing for the first time a high precision comparison between theory and data. The actual reach of 
the theory is affected by theoretical uncertainties associated to not having included higher order terms 
in perturbation theory, for which we provide an estimate, and by potentially overfitting the data, which 
we also try to address. Since in the EFTofLSS the coupling constants associated with the counterterms 
are unknown functions of time, we show how a simple parametrization gives a sensible description of 
their time-dependence. Overall, the fc-reach of the EFTofLSS is much larger than previous analytical 
techniques, showing that the amount of cosmological information amenable to high-precision analytical 
control might be much larger than previously believed. 
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1 Introduction 

The Effective Field Theory of Large Scale Structures (EFTofLSS) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 
14, 15, 16, 17, 18, 19, 20] provides the analytical framework that allows one to compute the distribution 
of dark matter and galaxies at large distances as a perturbative expansion in powers of the overdensity. 
So far, the EFTofLSS has been compared to simulation data for the case of the dark matter density 
power spectrum [2, 4, 9] and bispectrum [10, 11], the dark matter momentum power spectrum [9], the 
dark-matter vorticity slope [4, 21], the baryon power spectrum [4], the halo power spectrum and bispectra 
(including all cross correlations with the dark matter held) [13, 17], and the dark matter power spectrum 
in redshift space [22], The results have been very encouraging, showing that the EFTofLSS has a percent 
level agreement with the numerical data to a much greater wavenumber that formerly available analytic 
techniques (these analytic techniques are indeed incorrect if the EFTofLSS is correct). Maybe the most 
amazing result was that the EFTofLSS seemed to agree within roughly 2% with power spectrum data from 
the Coyote emulator [23], potentially up to the relatively high wavenumber of k ~ 0.6 /iMpc -1 [4, 9]. Very 
recently, while this paper was in advanced development, Ref. [20] appeared that attempted to analyze 
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with great precision the behavior of the dark matter displacement field in the EFTofLSS, using novel 
techniques that go beyond the simple power spectrum analysis. They find the EFTofLSS to fail against 
simulations at 1% level at k ~ 0.2/iMpc -1 , due to the appearance of stochastic contributions, with the 
EFTofLSS performing much better than former techniques. 

The EFTofLSS differs from former analytical techniques for two different reasons. First, the IR- 
resummation of infrared modes which was observed to be necessary to be treated non-perturbatively in 
order to have a well-defined perturbative expansion (see for example [24, 25, 26, 27, 28, 29]), is done 
in a radically different way than in these techniques, such as RPT [30]. According to a theorem by 
Frieman and Scoccimarro [31], generalized to the general relativistic context in [3, 32], the resummation 
of IR modes in the dark matter power spectrum should only affect the perturbative reproduction of 
the BAO peak, which appears in the power spectrum as oscillations in fe-space. Therefore, according 
to general relativity, a correct IR-resummation of this quantity should remove the residual oscillations 
in k-space between theory and data, without changing the UV reach of the theory compared to when the 
IR-resummation is not performed. This is achieved by the IR-resummation developed in the context 
of the EFTofLSS in [9], but, to our knowledge, not so in all formerly available techniques (including 
RPT [30] ). It should be stressed that the IR-resummation developed in [9] differs from former approaches 
in the actual implementation, which now respects general relativity, not in the conceptual fact that IR 
modes should be resummed to correctly reproduce the BAO peak, which had already been emphasized 
in different contexts [24, 25, 26, 27, 28, 29]. 

The second and most important difference between the EFTofLSS and other perturbative approaches 
is in the way short distance nonlinearities are treated. These other approaches, including RPT [30] or 
RegPT [27], assume (as SPT does) that the short-distance modes have a vanishing stress tensor. However, 
this is not an innocuous assumption: it implies that short-distance physics affects long distance dynamics 
only through the effect of the loops originating from perturbatively solving the nonlinear fluid equations. 
Rather, these loops receive a non-negligible contribution from modes so short that they are not in the 
perturbative regime. Even though short modes are not under perturbative control, they do affect long 
distance physics, and therefore need to be correctly parametrized. 

The EFTofLSS generalizes SPT by allowing for the most generic contribution of short modes at long 
distances. This results in extending the SPT equations to fluid-like equations, where the effect of short- 
distance modes at long distances is encoded in an infinite series of stress-tensor-like terms. The number 
of terms is infinite because all possible terms allowed by general relativity are introduced, in all powers of 
the long wavelength fields and number of derivatives. These terms are stress-tensor-like because they do 
not take the form that we normally have in a Navier-Stokes fluid, because the fluctuation fields include 
the tidal tensor of gravity, normally absent for fluids, and most importantly because the stress tensor 
depends on these fields in a manner which is local in space but non local in time [4, 6]. 

The expression of the effective stress tensor in a perturbative series of long-wavelength fluctuations 
is what makes the EFTofLSS the correct theory of the long-distance universe, in this superseding SPT. 
It however comes at a cost. Contrary to SPT, the EFT has in principle an infinite series of unknown 
parameters. However, the situation is not as tragic as it might appear at first glance. First, each of 
the terms of the effective stress tensor contribute only starting at a given order in perturbation theory, 
so that, to make finite-order predictions, only a finite number of terms are needed. In practice, all 
results obtained so far for dark matter have been obtained by using only one or two of these unknown 


3 


parameters 1 . Second, the parameters need to be measured in one of the following two ways. Either one 
can measure them directly by matching the predictions of the theory to long-wavelength observations (or 
to simulations, which are nothing but numerical experiments). This does not mean that the theory loses 
all predicting power, because each of the unknown coefficients comes with a specific functional form in 
wavenumber-space, so that not all information is lost. Indeed, this is the way we measure the Newton 
constant in general relativity 2 , or the way we measure F 7 r in the Chiral Lagrangian that describes pion 
interactions. The second way in which the parameters of the EFTofLSS can be measured is by measuring 
the effective stress tensor directly from dark matter particles, which, in contrast to the fluid elements, 
represent the correct degrees of freedom at short distances. This measurement can be done with small 
simulations, which only have to reproduce the nonlinear scale, and are therefore very fast and potentially 
more accurate. This method of measuring the parameters of the EFTofLSS leaves no free parameter when 
the theory is compared to data, and was pioneered in [2]. 

The reach of the EFTofLSS depends on the precision of the measurement of the paramaters, and 
therefore on the quality of the data available. Previous results were obtained from emulators such as those 
provided by CAMB [33] or Coyote [23], which have at least one-percent error, or with high cosmic variance 
simulations. These results seemed to show that one could could fit the data up to k — 0.6 h Mpc -1 staying 
within the error bars, but with an uncertainty that could push back the /c-reach as low as k ~ 0.4 h Mpc -1 
(see e.g. [4] or [9]). Furthermore, the objective of the first papers on the EFTofLSS was to focus on 
understanding the theoretical aspects of the theory, rather then to establish the precise fc-reach of a given 
fixed order calculation. 

In this paper, we compare the predictions of the EFTofLSS with data from a simulation (Dark Sky [34] , 
described in Sec. 2.1) with very small cosmic variance down to relatively low wavenumbers. Furthermore, 
since we have access to more detailed, quasi-direct, information about the data, it is possible to control 
many of the systematic errors that can occur in the comparison between simulations and theory 3 . We use 
a fitting procedure that properly incorporates the cosmic variance uncertainty on the power spectrum, 

1 Very explicitly, the prediction that agrees for the power spectrum at roughly percent level, at redshift zero up to 
k ~ 0.25 h Mpc^ 1 [2] at one loop, and up to k ~ 0.6 h Mpc -1 [4, 9] at two loops, and for the matter bispectrum at 
redshift zero at one loop up to k ~ 0.3 h Mpc -1 [10] is obtained using only one and the same parameter, c s (z = 0). 
Similarly, the prediction of the momentum power spectrum at one loop at redshift zero up to k ~ 0.3/i Mpc -1 [9] 
is done in principle by not only using c s (z = 0), but also c s (z = 0). In practice, c s {z = 0) can be inferred with the 
required precision at z = 0 by using an approximate scaling symmetry of the universe, so that c s (z = 0) might not 
be considered as a free parameter. Even if one were to consider c s (z = 0) as an additional parameter, one should 
consider that with this additional parameter the EFTofLSS is able to fit the dark matter power spectrum at all 
redshifts [16], and the same is expected to hold (thought it has not been verified yet) for the bispectrum and the 
momentum power spectrum. The prediction of the slope of the vorticity field does not require any new parameter. 
To similar precision, the prediction of the baryon power spectrum up to k ~ 0.6 h Mpc -1 at z = 0 requires one 
additional parameter [14]. 

2 General Relativity is indeed nothing but the Effective Field Theory of a massless spin-2 particle, and has 
therefore an infinite number of parameters. Luckily, and at the same time unfortunately, it is very hard to measure 
the effect of the additional terms. 

3 One of the authors would like to stress that it is a well known fact among people dealing with numerical 
(and experimental) data that to perform comparisons with exquisite precision, it is necessary to have comparable 
exquisite knowledge of the source of the data themselves. When previously existing data are subjected to new 
analyses, it can sometimes happen that systematic effects in the data that were previously irrelevant will become 
relevant to the new analyses, and we expect that data from TV-body simulations are no exception. 
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without trying to account for unknown systematics, finding that the prediction of the EFT at two loops, 
including only the lowest-order counterterm (associated with a single free parameter), agrees with the 
nonlinear measurements at redshift 2 = 0 up to k ~ 0.15 h Mpc -1 to within ^0.3% (the cosmic variance 
errorbar at that wavenumber). This represents a reduction in the fc-reach of the EFT with respect to 
former results, albeit with a much higher level of precision and accuracy of the numerical data. Note, 
in fact, that on the same numerical data, linear theory and two-loop SPT display a similar reduction 
in their k- reach, failing at k ~ 0.04 h Mpc -1 already by 1%, against the common knowledge that states 
that linear theory fails at 0.1 /iMpc -1 . Indeed, calculations in the EFT represent an expansion roughly 
in powers of /c/fcNL, so that a calculation at a given order will fail at lower and lower wavenumber as we 
increase the precision of the data more and more. Furthermore, as we describe in detail in the bulk of 
the paper, this reduction in the fc-reach is mainly due to a shift on the value of the speed of sound c^s, 
induced in turn by the more precise numerical data. As clearly highlighted in [16], the fc-reach of the 
theory obtained in former studies was indeed highly sensitive to the numerical value of c 2 s n y 

In previous work, we justified neglecting all other counterterms in the two-loop EFT prediction on 
the basis of order-of-magnitude estimates of their sizes. In this paper, we introduce a more detailed 
method for determining the importance of these counterterms, based on the fact that these terms should 
compensate for the contribution in loop corrections by modes with wavenumbers larger than the nonlinear 
scale. We find that the counterterms that have previously been included are not sufficient to make 
the calculation UV-insensitive (that is, to correct for the bulk of the UV contribution from the loop 
integrals), and that two more counterterms are necessary and sufficient for this purpose. With a total of 
three parameters, the EFTofLSS at two loops becomes UV-insensitive, and the A:-reach of the theory is 
increased to k ~ 0.34 h Mpc -1 , where the cosmic variance of the simulation is about 10~ 3 . 

We then investigate the impact of various other sets of counterterms, including a stochastic term, on 
the power spectrum prediction and its performance with respect to simulation data. We also extend the 
analysis to higher redshifts, finding similar results for the fc-reach, and study the time-dependence of the 
counterterms. Given that the number of modes scales as the maximum wavenumber cubed, the gain from 
using the EFTofLSS with respect to linear theory is still very large, somewhere between two and three 
orders of magnitude before considering the loss of information due to the marginalization over the three 
counterterms. Such an increase in the number of available modes might have huge consequences for our 
capabilities to explore the early universe in the next decade. 

2 Simulations and Fitting Procedure 

2.1 The Dark Sky simulations 

In this paper, we will compare the EFTofLSS prediction for the matter power spectrum to one of the Dark 
Sky series of dark matter-only N-body simulations 4 . These simulations use a Planck-like cosmology with 
cosmological parameters {D m , Ob, Oa, h, n s , as} = {0.295,0.0468, 0.705,0.688,0.9676,0.835}. Their initial 
conditions were set up using 2LPT, with an initial power spectrum generated by the CLASS code [35]. 

In particular, we will utilize the dsl4_a run of Dark Sky, which evolved 10240 3 particles in a box of side 
length 8/i _1 Gpc from an initial redshift z\„it = 93 to z = 0. The matter power spectrum was measured 
from dsl4_a snapshots at various redshifts on a Fourier grid with 8192 3 cells, with a constant shot noise 

4 http: / / darksky.slac.stanford.edu 
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contribution subtracted analytically, and averaged into bins of width A k = 2ir/L^ ox ~ 8 x 10 4 /iMpc 1 
The cosmic variance errorbar on the value of each bin is estimated by the standard Gaussian formula, 



where Ni is the number of modes contained in bin i. Where this number falls below 10~ 3 , we instead use 
5Pi/Pi = 10~ 3 as the errorbar, because this is roughly the precision of the theoretical computations that 
we will compare to the data. 

The combination of high resolution and large box size suppresses the cosmic variance of the power 
spectrum to sub-percent levels for k > 0.05/iMpc -1 , allowing a detailed comparison between theory and 
nonlinear measurements at much smaller wavenumbers than previously possible. Note, however, that 
this cosmic variance is already present in the (linear) initial conditions, rather than being significantly 
altered by nonlinear evolution. Therefore, when plotting various power spectrum predictions we replace 
the theoretical, averaged linear spectrum (generated by CLASS) by the spectrum measured from the 
initial conditions of dsl4_a (rescaled to the appropriate redshift). We make this replacement only for 
k < 0.1/iMpc _1 . By doing so, we can dramatically reduce the low-A; variance between the predictions 
and numerical data, as seen in the ratios we will plot later in this paper. 

We only make this replacement when generating the plots; the actual fits to the data (described in 
the next section) are performed using the theoretical linear spectrum. This is because it is more difficult 
to estimate the residual errorbars on the simulation measurements due to cosmic variance if the measured 
Pn is used in the theory prediction (or used to compute the loop corrections, which is also possible in 
principle), and the cosmic variance of dsl4_a is small enough that we can obtain satisfactory results 
without this extra complication (it would be interesting to pursue this approach in future work). 

2.2 Fitting procedure 

The EFTofLSS power spectrum consists of the usual SPT expansion plus a set of counterterms coming 
from the expansion and evaluation in perturbation theory of the dark matter effective stress tensor 
(on large scales). The parameter multiplying each counterterm is unknown and is fixed by the fitting 
procedure, but there are two complications that arise at this point. First, there is some ambiguity about 
which counterterms to include in a given prediction, and this ambiguity is only partially reduced by 
attempting to estimate the relative size of each term, since these estimates are only accurate at the order- 
of-magnitude level. In Sec. 3.3, we will describe a method to approximately estimate which counterterms 
should be included in a given calculation, while Sec. 4 will be devoted to an exploration of various 
combinations of counterterms that can be included, such as a stochastic term oc fc 4 . 

The second complication is that, apart from the aforementioned order-of-magnitude estimates, we 
do not know the maximum reach of the theory a priori. However, we can use the following fact to 
our advantage: as one increases k max (the maximum wavenumber used for a fit), we observe that the 
parameters change by remaining within the error bars of the fit (which shrink as we move to higher 
fcmaxi because as we increase k max we use more data points), and beyond some k max , the parameters start 
changing beyond the amount expected from the error bars. We interpret this as the fact that, by adding 
data, the parameters are better constrained and converge to their actual values, staying approximately 
constant as a function of k max . When k max enters the region where higher order terms become relevant, the 
best fit parameters change in a statistically unexpected way because the fitting method tries to compensate 


6 


Notation 


l max 

/'(it 


k 


reach 


Definition 

Maximum k of power spectrum measurements used in a fit 
For a given prediction, the maximum value of /c max for which the 
fit parameters are “stable” (as defined in the main text) 

The wavenumber at which the prediction fails, 
when parameters are fit using the region k < k^ t 


Table 1: Descriptions of various wavenumbers defined in Sec. 2.2 and used throughout the text. 


for the omitted terms. We therefore define k&t as the maximum h max for which the parameters have stable 
values and &: reac h to be the corresponding reach of the theory for these values of the parameters, that is, 
the maximum k for which the theory prediction is within the error bars of the data. 5 Note that fc reac h is 
not necessarily equal to k^. In principle, the theory curve could continue to match the data beyond k^ t . 
In practice, this does not happen because, by construction, fc reac h is the maximum h max beyond which 
the values of the parameters that would match the data change significantly. For the convenience of the 
reader, we summarize the definitions of h max , kfn, and fc reac h in Table 1. 

The stable region of the parameters can be identified by determining the 2cr region associated to the 
fit of each parameter at a given &: max . As we increase fc max , we expect the determination of the same 
parameter at the higher fc max to lie within the 2cr region of the lower h max , as otherwise the fit obtained 
with the new parameter up to the smaller /c max would be significantly worse 6 . We present figures of the 
parameters as a function of k mSLX later in the text, and we use this as a criterion to determine We also 
check that the stable region is always within a range of wavenumbers where the p -value of the comparison 
between theory and data is close to one. We believe that this method should help prevent overfitting. 
In the next section, each plotted power spectra corresponds to the best fit value of the parameters for 
the fcfit determined by this procedure, even though we could practically chose any feg t within the stable 
region, as the theory curve would not change much if the parameters were taken with a fcgt everywhere 
in the stable region. 

We estimate the theory error by plotting the values obtained by choosing one of the parameters to be 
la away from the best fit value that we have at 0.75fcfit, an d then re-fitting for the other free parameters in 
the prediction. Since physical results should be independent of the renormalization scale, this procedure 
encapsulates the effect of higher-order terms that are not included. This procedure is also affected by the 
smaller amount of data that we have at a lower k& t, which affects the determination of the parameters. 
Unfortunately, in an EFT, estimates of any results, such as the h-reach or the theoretical error, are most 

5 Since our goal in this paper is to perform a controlled comparison between the EFTofLSS and the Dark Sky 
simulations, we find it natural to use the uncertainty on the simulation output to define when the theory “fails.” 
Instead, if one is comparing to another simulation or to observational data, a different definition of “failure” could 
certainly be more appropriate. Note, however, that near-term observations are unlikely to require a precision as 
stringent as we require in this work, so in this light, our assessment of the failure of the power spectrum prediction 
can be seen as somewhat pessimistic. 

6 The criterion for allowing for a 2a discrepancy is a bit arbitrary. We could have chosen, for example, a ler 
discrepancy, even though we would find a bit odd to forbid fluctuations in the parameters of more than ler. Our 
estimate of the theoretical error will nevertheless be such that it will encompass the case if we had chosen a ler 
threshold. 


7 





reliably performed using the data themselves. One could use purely theoretical estimates such as those 
provided in [16], but they agree with the ones used here, given the fact that the theoretical error should 
be taken at the order of magnitude level. 

For numerical computations of the loop corrections in the EFTofLSS, we use a modified version of 
the Copter code [36] that makes use of the IR-safe integrands of [3] and the Monte Carlo integration 
capabilities of the CUBA library [37]. We use the same input linear power spectrum as the Dark Sky 
simulations, and require a relative numerical precision of 10~ 3 on all results. The IR-resummation method 
of [9] is then implemented in Mathematica, and is used for all EFTofLSS predictions presented in this 
paper. 


3 The Two-Loop Power Spectrum 

In this section, we discuss the two-loop matter power spectrum in the EFTofLSS. We first review the 
formulas that have previously appeared in the literature, describe procedures for determining the free 
parameters c 2 ^ and c 2 ^ associated with them, and present a comparison to redshift-zero measurements 
from the Dark Sky simulation. We then move on to re-examine the choice of which counterterms to 
include in the calculation, and find that considerations of UV-sensitivity necessitate the inclusion of two 
additional counterterms. The final subsection compares this three-counterterm prediction to data. 

In the conclusions of this paper, we will comment on possible numerical problems that might affect 
the precise comparison of the EFTofLSS with simulations. For the moment, we will proceed without 
questioning the reliability of the comparison we are performing. 


3.1 Formulas and renormalization 


We will begin by reviewing the one- and two-loop EFTofLSS predictions for the power spectrum as they 
have appeared previously in the literature. The one-loop EFT formula for the power spectrum is given 
by 

P E FT-i-ioop(fc, Z) = [Di(z)} 2 Pu{k) + [Di(z)] A Pi_\ oop (k) + P<g>(k, z) , (2) 

where 

p£!(k,z) = -2(27r) c 2 (1} ( z )[D 1 (z)] 2 ^P 11 (fe) , (3) 

while the two-loop formula is 


EeFT-2-1oop(^> Z ) = ^EFT-l-loop (k,z) + [Di(z)] 6 P 2 -loop(^) — 2(27 t)c^ 2 ) (z) - — -^P\i{k) 

«NL 

+ (2^)cl m (z)[D 1 (z)} 4 P^ op (k) + (27T) 2 ^1 + . 


(4) 


Expressions for P\-\oop(k) and P2-\oop(k) are given in [3, 4]. Expressions for (fc) are given in [4] for 

the £ = 2 case; the extension to C, ^ 2 is straightforward using the recurrence relations from [10]. The 
derivation of the ( k/k^p) A P\\ term, along with a discussion of the ( parameter and how it is related to 
the time-dependence of the counterterms, can be found in [16]. 

Implementation of these formulas in a comparison against numerical data requires the determination 
of the parameters c 2 ^ and c 2 s(2 ) ■ The procedure for doing so was essentially laid down in [4], modulo a 
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Figure 1: Ratio of the predictions of the EFT at two and one loops. By choosing a value of k ren small 
enough, we are able to make this ratio very flat, demonstrating that the higher order terms contained in 
-Peft- 2 -Ioop are negligible below at the wavenumber at which we impose the renormalization condition. 


few changes that we will highlight. If we only wish to use the one-loop prediction, Eq. (2), we only need 
the value of which can be found from the procedure described in Sec. 2.2: choosing the best fit 
value of c 2 ,y up to the k max for which it is relatively constant and the p- value is acceptable. One can also 
determine directly from the two-loop prediction (4), after c 2 ^ has been fixed in a way that we will 
describe below. This is the main procedure we will use in this work, for two reasons. First, the two-loop 
prediction will match the data over a larger region than the one-loop prediction, enabling a more precise 
determination of c 2 s yy Second, there is a larger risk of overfitting, and therefore of obtaining an incorrect 
value of if only the one-loop prediction is used (we comment further on this in Sec. 3.2.) 

We now describe how to determine . As discussed in [4], when we evaluate the two-loop diagrams, 
we find that, when the wavenumbers running in the loops are taken to be much larger than the external 
wavenumber k , these loops produce a contribution that is functionally of the form k 2 P\i{k). This con¬ 
tribution is degenerate with the lowest-order counterterm Pj^J , associated with c 2 ^. The simplest way 
to handle this degeneracy is to choose c 2 ^ in order to cancel the part of the two-loop terms that scales 
like k 2 Pn. To do so, we can use the fact that the term k 2 P\\(k ) C P %loop is the one that dominates as 
k —> 0 '. Therefore, we can determine c 2 ^ by choosing a very low wavenumber k ren , at which the terms 
that decay faster than k 2 Pn(k ) have become negligible, and require that 

F > EFT-l-loop(^ l 'reni %) — F > EFT-2-loop(^reiD Z) ^ ^s(2) — ^s(2)(^ren, ^s(l)) • (b) 

Any sizeable non-logarithmic dependence of c^ 2 ) on k ren should be taken as an indication that the 
renormalization scale has not been taken infrared enough for all the terms in P 2 -loop other than k 2 P\\{k) 
to be negligible. If k Ten is taken sufficiently low that logarithmic running of c 2 , 2 ^ is not expected (based 
on the slope of the linear power spectrum), then c 2 ^ should be essentially independent of k Ten . 

The determination of c 2 ^ can be made easier by a manipulation of the two-loop integral. Since 
any term that is of the form of k 2 Pu(k ) is irrelevant to compute, as it amounts to a redefinition of 
c ^( 9 ), we define a “UV-improved” version of P 2 -ioop by subtracting from P 2 -ioop a term which is a good 

'Locality in space, as well as momentum and matter conservation, forbid any term that is generated by modes 
higher than k to decay more slowly than k 2 Pn(k) at low fc’s. 
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approximation of its k 2 P\\{k) component: 


- 2fc 2 Fi. 


P 

(*) / 

J\a 


d\i Ah 2 lim 


2k I P 11 (k) 


i T-jintegrand / 7 —* —* 

where P 51 B (k,qi,-qi 


{«,«}>**» (2-) 3 (2tt)3 

, 02 , — 52 ) is such that the usual SPT diagram i 3 5 i(A:) is equal to 

Ai(&) 


f d qi d q2 ^integrand/; — -» -> — \ 

= /, ( 2^)3 ( 2^)3 P 51 ( A :, 91,-71,®,-52) 


( 6 ) 

(7) 


Using < V o ; mpr ° vcd) has the advantage that, if &: ren and fc m ; n are taken low enough, the that will be 
obtained from Eq. (5) will be very small, and in fact vanishingly small apart from running effects 8 . 

From this discussion, it is clear that, once c 2 ^ is determined from (5), the ratio Peft- 2 -Ioop/Feft-i-Ioop 
remains close to one up to a higher k than if a different value of c 2 ^ was to be chosen, as this ratio goes 

as 1 + EjdoopV-PEFT-i-ioop, which is much closer to one than 1 + ( k/k^ L ) 2 Pn(k)/PEFT-i-\oop ■ Fig. 1 
presents the ratio Peft-2-1oop/-Peft-i-1oop for two different values of c 2 ^ obtained by applying Eq. (5) at 
two different renormalization scales. The ratio has the expected behavior: it is quite close to one up to 
k ~ 0.1 /iMpc” 1 , where P^p causes it to deviate from one by ~ 1%. Indeed, this is the scale at which 

the one-loop EFT fails at one-percent level, as it lacks -F^Uoop^ furthermore, we have checked that as 
soon as we relevantly change the value of c^ 2 ^, the ratio deviates from one at a much lower k. Finally, 
we find that the ratio shows the same behavior whether we use the UV-improved version of T^-ioop or 
not; in the former case, we find c 2 ^ of order 0.02 (£:nl/(2 h Mpc ~ 4 )) 2 , while in the latter case we find 

c 2 ^.j ~ —2 (/cnl/(2 h Mpc ~ 4 )) 2 . The difference between these two values is the contribution to c 2 ^ that 
is removed by the UV-improved procedure. For the rest of this work, we will use the UV-improved version 
of P 2 -I 00 P) -^ 2 -Uoop mPr ° Ved ^, aR d choose c 2 s(0) = 0. This procedure can be thought of as choosing k Ten — k m [ n , 
and we will take k m i n ~ 5 x 10 -4 hMpc^ 1 , which is much lower than the lowest wavenumber we will use 
in our comparison between theory and data. 


3.2 Comparison to data at 2 = 0 

Before re-examining the self-consistency of the two-loop prediction in Eq. (4), we will first compare it to 
the matter power spectrum measured from the Dark Sky simulation at redshift z = 0, in order to enable 
an easier comparison with previous results. From Fig. 2, we can see that the UV reach increases as each 
set of higher-order terms is added: the linear prediction fails at fc reac h ~ 0.035 ZiMpc ^ 1 by ~ 1.4%, at 
one loop at /c rea ch — 0.08/iMpc^ 1 by ~ 0.5%, and at two loops at fc reac h — 0.15hMpc -1 by ~ 0.4%. In 
contrast, we see that there is no relevant improvement between linear theory and two-loop SPT (and all 

8 F 2 ( I V o r PrOVed) a ^ so ^ as a significant advantage in terms of computational cost, because the UV-sensitive term 
k 2 Pn(k) that is removed from P 2 _i 0 op is typically much larger than the UV-insensitive term, P^Uoop^ we are 
interested in calculating, in the range of wavenumbers relevant for this work. Even in -P 2 -Y 00 ]^ mproved ) j there will 
be some residual UV-sensitive terms, but these are also much smaller than P^Uoop^ Therefore, for a given level 
of accuracy that we desire for -F^Uoop^ a smaller number of evaluations of the integrand (around one order of 

magnitude fewer, in fact) is required when computing ^ 2 ^ 00 ™^°™^■ The algebraic expression for P^-Voop 1 " 1 "° Ved ^ 
is much more complicated than for P 2 -ioo P , causing the UV-improved integrand to be ~10-20% slower to evaluate 
than the regular integrand, but this does not interfere with the significant computational gains of l ^"p mproved '* ■ 
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Figure 2: The two-loop EFTofLSS prediction for the z = 0 power spectrum, when one includes only 
one counterterm (associated with the speed of sound c^), along with various other theory predictions. 

The EFT curves use a value of c 2 ^ ~ 0.53 (&nl/(2 /tMpc -1 )) 2 . We can see that the theory performs 
better and better as higher order contributions are included. The blue shading represents the variation 
of the result if we perform the fit to determine c 2 ,^ up to 0.75/cfit and choose the two values lu away 
from the central value, where fcgt is the wavenumber beyond which c 2 ^ begins to deviate from the value 
determined at lower k. For k < 0.1 h Mpc -1 , the linear power spectrum in the theory prediction is 
replaced with the power spectrum measured from the initial conditions of the simulations, allowing for 
a dramatic reduction in the variance of the i'theory/-f > NL curves at these wavenumbers, but also implying 
that the cosmic variance errorbars (represented by the grey shading) do not reflect the uncertainty on 
the curves for k < 0.1 h Mpc -1 (as discussed in Sec. 2.1). The fc re ach of the EFT at two loops is about 
k Te ach — 0.15 h Mpc -1 , where the cosmic variance is about 0.4%, even though there is large theoretical 
uncertainty. The fc-reach is smaller than what was previously presented in [4], where the errorbars were 
taken to be ~2%, because the much higher precision of the available numerical data allows the choice of a 
lower k ren , which eliminates the strong cancellation between various two-loop terms that was seen in [4], 



analytic techniques prior to the EFTofLSS such as RPT and RegPT, which differ from SPT only by the 
resummation of the IR-modes, which are irrelevant for the UV reach of the theory 9 ). The blue shaded 
region in Fig. 2 represents the difference between fitting up to fcg t or instead fitting up to 0.75 and 
choosing the values of c 2 ,^ lcr away from the best fitting point; this represents a rough estimate of the 
theoretical error associated with the prediction of the EFTofLSS at this order, and should be taken at 
the order of magnitude level. 

We now make a few comments on these results. The first is on the importance of the IR-resununation. 
Without performing the IR-resununation, the EFT prediction would be off with respect to the data by 

9 This point is quite unappreciated in the literature, so we repeat it here. RPT is sometimes claimed to improve 
the UV reach of the theory. However, to our understanding, it is supposed to be just an IR-resummation, and 
therefore if it improves the broad-band UV reach of the theory it violates General Relativity. It is therefore 
incorrect and should not be considered as a way to increase the broad-band UV reach. Alternatively, one should 
consider RPT as a fitting function. See [14] for a more detailed discussion. It should not be forgotten that the 
fact that RPT violates General Relativity (and cannot therefore be a correctly implemented IR-resummation) was 
already pointed out in the original RPT paper [25], whose focus indeed was not on the broad-band fc-reach. 
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Figure 3: Same as Fig. 2, but using the procedure from previous papers (e.g. [4, 9, 16]) to fix c 2 ^ and 
c s( 2 ) ( see the ma i n text for details). If we allow for a uniform 2% error budget on Rnl ( as was required 
in previous work, due to the use of the Coyote emulator), and use a k ren comparable to previous work 
(fcren 0.23 /iMpc *), we find results that are consistent with [4, 9, 16] (the ~1.5% offset between the 
one- and two-loop curves was actually about 0.7% in the cosmologies and data considered in [4, 9, 16]). 
One sees that at low k, there is an offset between theory and data. Such an offset cannot be ruled out by 
using power spectra from the Coyote emulator, since its output has a systematic errorbar of at least 1%. 
Instead, for Dark Sky, we cannot allow for that level of systematic error, and forcing the absence of the 
offset at low k affects the k -reach of the theory, leading to the results in Fig. 2. 


oscillations of order ~ 2% (see [4]). Due to the smallness of the current error bars, it would therefore be 
impossible to impose the theory to have a good fit to the data, and so we see that IR-resummation is 
essential to performing this high-precision comparison. 

Next, we note that we have chosen a smaller k ie n than in previous works, starting from [4], allowing 
for a smaller theoretical error. Notice that since the renormalization procedure to determine c 2 ^ depends 
only on the theory calculation, it would have been possible to choose a low k ren also when comparing 
with more noisy data, such as the output of the Coyote emulator. This procedure is correct, but it has 
the technical inconvenience that, as we explain in the next paragraph, one would have not been able 
to compare the value of c 2 ^ at one loop and at two loops, which is an interesting, though delicate, 
consistency check of the theory, that was performed in [4]. This procedure is also more sensitive to the 
evaluation of the numerical loops and IR-resummation, as the overall effect of k 2 P\\ drops rapidly. 

Overall, it would be incorrect to interpret the higher fe-reach of [4] as simply due to a choice of too 
large a value of k ren . The renormalization scale chosen in [4] is actually not particularly high, when we 
consider that the low- 7 offset between Reft- l-loop and -Peft- 2 -Ioop at the k re n chosen in [4] is less than 
0.5%, which is not large given the error bars on the data in [4]. The reason why this choice leads to a 
sizable effect on the /c rea ch of the theory is due to the fact that, as explained in detail in [16], the value 
of c 2 ^ that happens to be chosen in this procedure leads to an accidental cancellation between R^Uoop’* 

and P l-irx) P • When combined with the poor precision of the numerical data used in [4] , the k ie ac h is made 
larger. For reference, in Fig. 3 we show what happens when we perform fits to the Dark Sky data using 
the procedure from previous papers [4, 9, 16], also allowing for a larger error budget on the nonlinear 
power spectrum. 

We also wish to highlight the importance of the two-loop contribution in determining from the 
power spectrum. Without the inclusion of the two-loop terms, one could choose the value of in 
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such a way as to potentially match the data up to k ~ 0.11 hMpc -1 . If c 2 ^y is instead determined from 
the full two-loop prediction, the value of we find, ~ 0.53 (&nl/(2 h Mpc -1 )) 2 , is roughly 30% 
smaller than that obtained from the one-loop fit, — 0.75 (&nl/(2 h Mpc -1 )) 2 , and causes the one- 
loop prediction to fail at k — 0.08 /iMpc^ 1 . This reveals that using the one-loop prediction on its own 
potentially leads to overfitting. 

Indeed, in the limit of infinitely precise data, we would like to determine c 2 ,y at as low a wavenumber 
as possible. In practice, though, there is a minimum wavenumber below which the effect of the leading 
counterterm becomes smaller than the uncertainty on the data, and this wavenumber will act as a lower 
bound for the region where it safe to fit for c 2 ^y This minimum wavenumber is around k ~ 0.06 h Mpc -1 
(estimated by when P£ T s ee /Pn equals the uncertainty on our data). Meanwhile, the one-loop prediction 
begins to fail around k ~ 0.08 h Mpc -1 , but fitting up to this point would introduce a substantial bias 
on the resulting value of c^, since, by our definition of “failure,” at k ~ 0.08 /i Mpc^ 1 the two-loop 
terms contribute by an amount equal to the uncertainty on the data, which in turn is similar to the 
size of the effect of at these wavenumbers. Therefore, the region in which it is safe to fit for c 2 ^y 
using -Peft- 1 -Ioop alone cannot be larger than roughly 0.06 h Mpc -1 < k < 0.07 /iMpc -1 , but it would be 
impractical to try to fit over such a small region. On the other hand, our procedure of fitting for 
directly using the two-loop prediction allows us to use a much larger range of wavenumbers. 

Furthermore, as mentioned earlier in this section, another contribution to the mismatch of c 2 ^ from 
the best fit we obtain at one loop and at two loops could be associated to the fact that our two-loop 
renormalization procedure can be thought of as determining c^ 2 ) at low wavenumbers, and at around 
k ren ~ breach at two loops. The one-loop renormalization procedure instead determines around the 
breach °f the one-loop computation. If the two-loop contribution were to have a logarithmic running at high 
wavenumbers, the value of that we obtain by renormalizing the two-loop theory at high wavenumbers 
would be different than the one obtained by renormalizing at low wavenumbers, as it would need to absorb 
the running of c^ 2 ) between the two renormalization scales, which we do not account for. There are also 
effects due to the fact that the measurement of is done at a non-vanishingly small value of k ren , and 
that vanish only for k ren —> 0. Overall, this mismatch has no effect on the two-loop result present in this 
paper (because it is degenerate with the additional counterterms), nor on the one-loop result as well if 
we allow for to have different values at one and two loops. Since, due to the non-scale-free nature of 
our universe, precisely determining if the one-loop or two-loop diagrams have a logarithmic contribution 
is very hard if not impossible, and given the fact that this discussion is relevant only for comparing one- 
and two-loop results, not for the maximal reach of one- and of two-loop calculations, we do not explore 
this issue further. 


3.3 UV-sensitivity of the standard two-loop calculation 

The results of the previous section show that in order for the theory to agree with data beyond k — 
0.15 hMpc -1 , we need to add more counterterms or loops to the two-loop formula given in (4). Following 
the notation of [10], the next terms that are supposed to contribute to the effective stress tensor ( dr) Pl l 
are 

• terms quadratic in the fields, which (as discussed in App. A) take the form 

(dr) pi l D (1 - 5) x {<T ( -n(l)f - $) ■ 91 [d 2 ^] 2 • d l [&> & $ djd k (j^ , d'd^djd 2 ^ ; (8) 
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• higher derivative terms: 


(3rV 3 ; 


(9) 


• stochastic terms: 


(. dr) Pl i d PAt . 


( 10 ) 


(See [2] for the definition of At.) 
• cubic counterterms: 


(dr),; D #5 3 ,.... (11) 

Since there are many of those, we just wrote a representative one. 

Recall that the two-loop formula in (4) was derived from a form of the effective stress tensor that 
neglects all of the above terms, ft might seem particularly strange for a two-loop calculation not to include 
the quadratic terms. The justification provided in [4] was that the quadratic counterterms are included 
to remove the UV-sensitivity of diagrams that are not enhanced by as many factors of (2n) as and 
c s( 2 ) are - This reasoning is only partially correct, as shown in [10] in the context of the bispectrum. If 
we focus on the one-loop bispectrum for simplicity, at the order at which we are working, only three of 
the terms in (8) give rise to non-degenerate bispectrum shapes. In this case, if we take for simplicity a 
no-scale universe with slope n = — 1, where the corresponding diagrams are logarithmically divergent, 
we find that the resulting three counterterms remove divergences that are all enhanced by one factor of 
(2vr), but, for two of of the three coefficients, they are suppressed by numerical factors that undo the 
enhancement by (27r) 10 . Therefore, this type of 27r-counting argument turns out not to be particularly 
robust, and so we look for a more accurate strategy to assess which terms should be included at a given 
order in perturbation theory. Indeed, we can make a stronger distinction between which counterterms 
should be included in a two-loop calculation by examining the UV-sensitivity of the two-loop integrals 
directly, as we are going to explain next. 

Before describing the relevant procedure, let us enumerate the counterterms in the power spectrum 
that result from the stress tensor terms we listed above: 


• quadratic counterterms: 


Tquad. counterterms(^j *) 

(2t r 


( 12 ) 


&Nl/ 


(cow drrV)+drT’m+RC 3 ’®), 


where -Pm^’ 0 ’ 1,2,3 )(A;) can be found using the techniques and the kernels in the Appendix of [10], 
and we have kept only the contribution that is not degenerate with k 2 Pu(k ); 

10 The reason why this happens, also at finite momenta, is that an integrand of the form f d 3 q q^j is not ro- 
tationally invariant, but it can be written as ^ f d 3 q q 2 . So, even though there is a factor of (27 t), the resulting 
enhancement is reduced by the numerical factor (1/3) originating from the non-rotational-invariance of the inte¬ 
grand. In the case of the three quadratic counterterms, one diagram is numerically enhanced to reduce such a 
suppression. 
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• higher-derivative counterterm: 


Pl-deriv. counterterm(^i %) — 2(27t) D\(z) £4(2) 


—V 

^NL / 


Pu(k) ■ 


(13) 


• stochastic counterterm: 


Pstoch(^) z) 


(2tt) 2 D^z) 2 c stoch (z) 


k \ 4 1 

&nl ) A:nl 3 


(14) 


cubic counterterms: they are degenerate with k 2 Pn(k), and therefore do not need to be included 
in the calculation. 


The time-dependence of each term is accounted for by the growth factor D\(z) and by the explicit time- 
dependence of the coefficients. The factors of (2n) have been chosen in such a way that these terms have 
the same numbers of (27r)’s associated with the two-loop diagrams that generate these terms. 

As previously stated, the role of the counterterms in the EFTofLSS is to make the result of a calculation 
insensitive to the loop integrals that are performed in perturbation theory when the integrands are 
evaluated at high wavenumber where perturbation theory does not apply, and to instead parametrize the 
correct contribution of short distance physics at large distances. This suggests that a way to check that 
our two-loop calculation is consistent is to ask how much it depends on the region of the two-loop integral 
evaluated on momenta larger than the nonlinear scale. If the result is UV-insensitive, we should obtain the 
same result at low k if we evaluate the loop integrals with cutoff A = 0.68 hMpc -1 (taken to be a rough 
proxy for the nonlinear scale at z = 0, which we estimate to be about 2 k reac h(z = 0) ~ 0.68 h Mpc -1 ) or 
A = oo 11 by simply readjusting c^. Equivalently, if we use the UV-improved integrand, we should find 
the same result. 

Upon numerically evaluating P 2 -ioop at these two cutoffs, we find that, to the relevant level of precision 
required by the simulation data, the result is indeed not the same, and that additional counterterms must 
necessarily be added to make the result UV-insensitive to the precision we require. In particular, we can 
obtain an estimate of the size of the additional counterterms that are required in order to make the result 
UV-insensitive by fitting the two-loop result integrated from A = 0.68 hMpc -1 to 00 with the additional 
counterterms mentioned above. In formulas, we impose 


P, 


(UV-improved) 


2-loop 


A=0.68 h Mpc 


- p 


(UV-improved) 


2-loop 


A=c 


(15) 


£4 4 DV> AC' ‘V) + 2 ( ^4 UV) {£) «,(*) + WC (£)^ 


We add only P^lt' ° and not any of the other three quadratic counterterms from (12) because the shapes 


1-loop 

of these terms are quite similar. 

We present in Fig. 4 the ratio of 


dition or not of two sets of counterterms: {c^ 


over P 2 _ lo 
A=oo ^ 

(UV) (UV) „(UV) 

> c 4 


(UV-improved) 


A=0.68 h Mpc 


, with the ad- 


c stoch} an d {cS UV) ,4 UV) }- We can see that 


without the addition of the new counterterms, P^-Yoop™^ 0 ™^ 


A=0.68 h Mpc 


is significantly different from 


In practice, we use A = 60 h Mpc 1 for the numerical evaluation of the integral. 


15 












(UV-improved, A=oo) . , (UV-improved,A=0.6B * Mpc 
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Figure 4: The red dashed line shows the ratio of the calculations of i^-ioop with cutoff A = oo and 
A = 2 fc reac h(z = 0) = 0.68 /iMpc' 1 , while the blue dotted and black solid lines show the same ratio 
but adding the two or three counterterms to the A = oo FV loop calculation, respectively. We see that 
the difference between the A = 0.68 /i Mpc ^ 1 and A = 00 calculation of P 2 -I 00 P can be absorbed by the 
counterterms. This is an important consistency check of the EFTofLSS, indicating that the unknown 
short-distance physics affecting the loop corrections can be accounted for by the EFT counterterms. It 
also tells us that the two-loop power spectrum in the EFTofLSS should minimally include the counterterms 
corresponding to c\ and C 4 in order to be insensitive to our assumptions about the UV behavior of the 
theory. 


Pr 


(UV-improved) 


2-loop 


A=c 


However, after the addition of the counterterms, the result is insensitive to the contri¬ 
bution of the integral from A = 0.68 /iMpc ^ 1 to 00 , indicating that a sensible two-loop computation must 
minimally include the c\ and C 4 counterterms. (Other combinations of the cp, C 4 , and c s t oc h counterterms, 
not shown, fail to absorb the UV-sensitivity to an adequate level.) The fact that that the counterterms of 
the EFTofLSS are able to absorb this UV-sensitivity is a nice demonstration of the internal consistency 
of the theory 12 . 

The best-fit numerical values for the coefficients c ( | UV \ C 4 UV \ and give us a sense of the numerical 
value of the coefficients that we expect to be generated by the uncontrolled UV physics. We find the 


12 From a field theoretical point of view, one can argue that the correctness of the EFTofLSS is manifest already 
from its construction. 
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following values when all three terms are included 13 

c ; uv) = - 1.1 (fc NL /(2 h Mpc -1 )) 2 , (18) 

4 UV) = -5.3 (/cnl/(2 h Mpc ” 1 )) 4 , 
c stoch = 6-2 x 10 4 (fc NL /(2/r Mpc -1 ))” , 

and the following values when the stochastic term is omitted: 

c ( uv > = —1.6 (fc NL /(2 hMpc -1 )) 2 , (19) 

4 UV) = -7.0 (&nl/ (2 /iMpc -1 )) 4 . 

The fact that the values of c^ UV ^ and are not strongly affected by the presence of the stochastic 

term is further justification that the majority of the UV-sensitivity of i^-ioop can be removed by these 
first two terms. Only these two counterterms need therefore to be necessarily included in order to make 
the two-loop calculation UV-insensitive. 

Fig. 5 shows the relationship between the terms we add to the effective stress tensor and a schematic 
representation of the lowest-order loop diagrams (in either the power spectrum or bispectrum) they are 
associated with via renormalization. Note that the fact that the c 2 s ^d 2 5 term renormalizes P 13 , one of 

the terms in Pi_i oop , implies that the C 4 < 9 4 h term renormalizes not just P 2 -loop, but also and Ppp^p; 

since, roughly speaking, these two terms behave similarly to Puioop but with an extra factor of k 2 . This 
implies a possible physical correlation between the values of C 4 and c\, and between C 4 and c 2 ^, but we 
will not explore this further in this work. 

3.4 The UV-insensitive prediction at z = 0 

As we have just argued, in order to make the two-loop calculation UV-insensitive, it is necessary to add 
two more counterterms: one of the quadratic terms, and the higher-derivative term. The self-consistent 
two-loop prediction therefore has an overall number of three parameters to be fixed from observation or 
simulations at a single redshift. 

13 If instead we take A = 1 x peach (2 = 0) we obtain 

4 UV) = -2.7(fc N L/(2/iMpc - 1 )) 2 , (16) 

c| UV) = -13 (pm/(2 h Mpc -1 )) 4 , 
c<£2 = 13 - 5 x 104 (^nl/(2 h Mpc -1 )) 7 , 

and the following values when the stochastic term is omitted: 

c[ UV) = -3.9 (fc NL /(2 h Mpc " 1 )) 2 , (17) 

ci UV) = -17(fc NL /(2hMpc- 1 )) 4 . 

These numerical values should be taken as a rough over-estimate of the size of the counterterms that is required 
to make the prediction UV insensitive. We clearly expect the strong coupling scale of the theory to be larger than 
1 x fc reac h- But, given our lack of knowledge of the precise value of this scale, we present values assuming that the 
strong coupling scale is either 1 x p ea ch or 2 x p eac h to give a rough interval for the expected numerical values of 
the counterterms, an interval that we hope is large enough to include the true values. 
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Figure 5: The correspondence between the terms we add to the effective stress tensor and a schematic 
representation of the lowest-order diagrams (in either the power spectrum or bispectrum) they are asso¬ 
ciated with via renormalization. 


In Fig. 6, we present a comparison of this prediction to simulation data at z = 0, using the techniques 
of Sec. 3.1 to fix c^ 2 n and applying the fitting procedure of Sec. 2.2 to the determination of c^, c\ and 
C 4 . We find that the /c-reach is brought up to k ~ 0.34 /iMpc -1 , where the cosmic variance is ~ 10 -3 , 
with a theoretical uncertainty that could bring the reach down to about k 0.26 hMpc 1 . The size of 
the cosmic variance shows the remarkable level of precision of the comparison. The parameters that we 
find are given by 

c 2 s{1) ~0.48(fe N L/(2hMpc- 1 )) 2 , (20) 

ci ~ -0.74 (&nl/(2/iMpc - 1 )) 2 , 
c 4 ~ -6.4 (fc NL /(2hMpc _1 )) 4 , 

and they are consistent with the sizes that are expected from the UV-dependence of the calculation. 

As anticipated in Sec. 2 . 2 , an important check that we carry out in order to ensure that we are not 
overfitting is to ensure that the parameters that we determine from the fit are constant as we increase the 
fitting region. We present the results in Fig. 7. In the upper figure, we can see that when we use all three 
counterterms, the parameter obtained from the fit is constant over the range k < 0.32 /iMpc -1 . 
Here by constant we mean that the curve is constant within the error bars (shaded) determined by the 
fitting procedures. We also present curves for the value of obtained when only a fraction of the 
necessary counterterms are included. In the blue curve, no additional counterterms are included, and the 
curve is never flat. When either ci or c 4 are included, the curve is flat until k — 0.22 /iMpc -1 , when 
then it starts deviating significantly from being flat. This can be understood by noticing that at low fe’s, 
the higher order counterterms are not very important, so that can be determined very well without 
the complete set of the relevant counterterms being included. However, at k > 0.22 /iMpc ^ 1 these terms 
start to become important, and it is impossible to have a flat curve for because the term with 
is trying to compensate for the lack of a relevant counterterm. The curves for the parameters c\ and c 4 
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Figure 6: The prediction of the EFTofLSS at two loops after the inclusion of both one of the 

quadratic counterterms, ~ k 2 Pi_\ oop , and the higher-derivative counterterm, ~ k 4 Pn(k). We see that 
the fc-reach is significantly improved (compared to the prediction including only one counterterm) to 
breach — 0.34/iMpc -1 , where the cosmic variance is about 10 -3 . The theoretical error, shaded in purple, 
is estimated to potentially decrease the fc-reach to fc rea ch — 0.26 h Mpc -1 . As discussed in Sec. 3.3, we 
consider this to be the most theoretically justified calculation of the EFT at two-loop order. 


follow a pattern similar to the one of c 2 ^, being flat until k — 0.32 /iMpc -1 . The only difference is that 
the error bars shrink relevantly for k > 0.22 h Mpc -1 because, as for the case of the c 2 ^ curve, these 
counterterms begin to be sizeable at these wavenumbers. We conclude that if we take k&t = 0.32 h Mpc -1 , 
we are probably not overfitting the data. More checks on not overfitting the data are presented in Sec. 5.2. 

4 Other combinations of counterterms 

The prediction we present in Sec. 3.4 contains the minimum number of free parameters required to render 
the calculation UV-insensitive at the level of precision we are concerned with in this paper. Nonetheless, we 
can also ask what happens when we remove or add various counterterms to this prediction. This is useful 
for gaining insight into the effect of different combinations of terms, and for performing consistency checks. 
It is also useful for discovering whether there is a version of the prediction with fewer free parameters 
that performs just as well; in fact, even if we can estimate that some counterterms are necessary to be 
included to make the calculation UV-insensitive, it could well be, at least in principle, that, once we send 
the cutoff of the calculation to infinity, the finite contribution of the given term happens to be small and 
need not be included, allowing for a prediction with fewer parameters. 

4.1 Quadratic counterterms: ~ k 2 Pi.\ oop 

We begin by adding the quadratic counterterms in (12) to the one-counterterm two-loop EFT prediction 
in (4). The functional form of each of these terms is quite similar, so we study the match of the theory to 
data first when one term at a time is added, then two or three terms. The results are presented in Fig. 8. 
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Figure 7: From top to bottom, values of the counterterms c^, c\ and c 4 at z = 0 as obtained from 
our fitting procedure as a function of the k max of the fit. The results are presented for various choices of 
the counterterms being included. In shading is the 2 a errorbar from the fitting procedure, with the la 
errorbar for the three-counterterm fit shown in long dashed lines. The presence of a flat region in k mSLX is 
interpreted as suggesting that a certain parameter is being well measured and the A; m a X of the fit has not 
been overestimated. When all counterterms are being used, we notice the presence of a flat region for all 
of the three parameters, ending at /c max ~ 0.32 /iMpc -1 . We conclude that the fcg t should be taken to 
be ~ 0.32 h Mpc -1 at z = 0. 
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- 2-loop EFT, c 2 (1) = 0.53 (A' N l/[2/j Mpc *]) 2 

— withP» c 2 (1) = 0.33 (W[2/J Mpc- 1 ]) 2 , c 0 =3.5(W[2/iMpc- 1 ]) 2 

— withP^op 1 ’- c 2 (1) = 0.48 (W[2/J Mpc- 1 ]) 2 , Cl = 0.47 (W[2/i Mpc" 1 ]) 2 

■ with P ( “, c 2 (1) = 0.5 (k NL /[2h Mpc" 1 ]) 2 , c 2 = 1.9 (k Nt J\2h Mpc" 1 ]) 2 

■ with «p 3> . c 2 (1) = 0.51 (k m J\2h Mpc- 1 ]) 2 , c 3 = 3.5 (k N[ J\2h Mpc" 1 ]) 2 


- with pXp 0 ’ + ^noop 1) . c ?(i) = °- 65 (*nl/[2 h Mpc" 1 ]) 2 , c 0 = -4.9 (* NL /[2/t Mpc" 1 ]) 2 , Cl = 1.3 (k N J[2h Mpc" 1 ]) 2 

with<op 1) + ■ p( iToop 2) . c 2 (1) = 0.61 (W[2/t Mpc" 1 ]) 2 , Cl = -2.4 (W[2/t Mpc" 1 ]) 2 , c 2 = 12 (W[2/i Mpc' 1 ]) 2 

- with pXp 2 ’ + ^noop 3) . c 2 (1) = 0.68 (k NL /[2h Mpc" 1 ]) 2 , c 2 = -16. Mpc" 1 ]) 2 , c 3 = 26 (*W[2A Mpc" 1 ]) 2 



k [h Mpc 1 ] 

Figure 8: Prediction of the EFT at two loops after the inclusion of the terms in (4) plus various combi¬ 
nations of the quadratic counterterms. We can see that if we add just one of the counterterms, the fc reac h 
is increased to fc reac h — 0.23 — 0.26 h Mpc -1 , depending on the counterterm that is chosen, where the 
cosmic variance is about ~ 2 x 10~ 3 . An estimate of the theoretical error is shaded in purple, showing 
the possibility of the decrease of the fc reac h all the way to fc reac h — 0.18 /iMpc^ 1 . When we include two 
of the quadratic counterterms, the reach of the theory can be further increased to /c rea ch = 0.3 h Mpc -1 , 
although this is likely the result of overfitting, as described in the main text. 


When we add just one of the counterterms, we see that the EFTofLSS matches the data up to 
breach — 0.23 — 0.26 /iMpc^ 1 , depending on the counterterm we use, where the cosmic variance of the 
data is as small as 2 x 10~ 3 . The sizes of the numerical prefactors when each term is included separately 
are given respectively by 

c 0 ~ 3.5 (fc N L/(2hMpc _1 )) 2 , (21) 

ci ~0.47(fc NL /(2/iMpc- 1 )) 2 , 
c 2 ~ 1.9 (fc N L/(2/iMpc _1 )) 2 , 
c 3 ~3.5 (fc N L/(2^Mpc" 1 )) 2 , 

while c^/jj is in the range 0.33 (&nl/(2 h Mpc -1 )) 2 < c^ 1( < 0.51 (&nl/(2 /iMpc' 1 )) 2 , with slight vari¬ 
ations as we add more terms. As expected from the analysis of the bispectrum [10], we see that for 
these terms to be relevant, their coefficients need to be boosted with respect to the value of c 2 ^ by large 
factors. The value of ci is therefore compatible with what is expected from the UV sensitivity of Th-ioop 
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in (18) 14 . 

When we add two of these quadratic counterterms at the same time, we see that the reach of the 
theory can be boosted to as much as k Tea c h = 0.3/iMpc -1 . However, we interpret this higher match of 
the theory with data as an overfit of the theory. In fact, in some cases, this higher reach is achieved by 
boosting the coefficients of the counterterms to what we interpret to be very large values, in such a way 
that the two contributions cancel each other, when on the contrary these contributions are not expected 
to be canceling against each other. This interpretation is confirmed by the fact the numerical coefficients 
that we obtain when we add two counterterms become quite large with respect to c 2 ^ (this is particularly 
true for C2 and C3) and to what is roughly expected from the UV in (18); for example, when ^ 

and -Puoop 2 ' ) are both included, we find that 

c 2 s{1) ~ 0.61 (A: NL /(2/iMp C - 1 )) 2 , (22) 

ci ~-2.4 (/c NL /(2/rMpc _1 )) 2 , 
c 2 ~ 12 (fc NL /(2/rMpc _1 )) 2 , 

and most importantly they change a lot with respect to their numerical values when only one term 
was used. Furthermore, they take opposite signs so that they can indeed cancel each other. In 
the case of Cq and Cj, it is mainly the mismatch from what is expected from the UV that pushes 
towards the interpretation of the increase of the fc reac h as an overfit. These considerations tell us 
that to ensure that we are not overfitting the data we should not rely only on the consistency of 
the measurement of the parameters as a function of but also on the estimated size of a term 
from the UV and on the change of the value of the terms as we include additional parameters. 

When we add three counterterms, and we impose that the value of c 2 ^ is not changed by 
more than a factor of two with respect to the value that we find when not including these terms 
(a fact that is justified by the hierarchy of the various contributions), we find that the h-reach 
of the theory is not relevantly improved, so we have not plotted these curves in Fig. 8. We 
therefore conclude that for the precision of the given data, and restricting only to the quadratic 
counterterms, it is enough to include only one of them, where the reach of the theory is boosted to 
about breach — 0.23 /rMpW 1 . This breach is inferior to what is obtained in the consistent calculation 
that was presented earlier, where we include also the higher derivative counterterm. 

Finally, we make an additional comment. As noticed in [4, 6], the EFTofLSS is local in space, 
but non-local in time. This means that the term Pi_^ op should actually correspond to three different 
terms, with slightly different functional forms (see [16] for the most clear presentation). We have 
checked that the functional form of each of these terms is however highly degenerate with the one 
from the quadratic counterterms and the k 2 Pu{k) term; and in fact, when we use them in the fit to 
replace the quadratic counterterms, the result is not significantly different. Similarly, in the non- 
local-in-time treatment, each of the quadratic counterterms leads to two independent functional 

14 As we stressed, the estimate of the induced size of a counterterm from UV physics should not be overinterpreted 
as more than what it is, a rough indication of the expected numerical value. In particular, we argued that C4 should 
also be included in the calculation, but its inclusion is expected just to provide an order one correction to the size 
of Ci obtained from the fit to the data. 
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- 2-loop EFT, c; (1) = 0.53 (k NL /[2h Mpc d) 2 



k [h Mpc 1 ] 

Figure 9: The prediction of the EFTofLSS at two loops after the inclusion of the higher-derivative 
counterterm k' i Pn(k). We see that the fc re ach is improved to fc reac h — 0.24 h Mpc -1 , quite similar to the 
increase we saw in Sec. 4.1 when adding a single quadratic counterterm. 


forms, but we have checked that they are also degenerate with the quadratic counterterms in the 
local-in-time approximation and the k 2 Pn(k ) term. 

4.2 Four-derivative counterterm: r^j k 4 Pn 

We next pass to study the effect of including only the higher-derivative counterterm (13) in the 
two-loop calculation with the c 2 ,^ counterterm only. The result is presented in Fig. 9. We see 
that when we add this term alone the fc-reach of the EFT is boosted to k ~ 0.24 h Mpc -1 . The 
parameter values that we find are 

c 2 s(1) ~ 0.48 (fc NL /(2/rMp C - 1 )) 2 , c 4 ~ -2.6 (fc NL /(2 h Mpc -1 )) 4 , (23) 

which are compatible with what is expected from (18). We see that the inclusion of this term 
results in a fc reac h that is very similar to that obtained from adding a single quadratic counterterm 
in Sec. 4.1. Evidently, either a single quadratic term or a single higher-derivative term are not 
sufficient to reproduce the results of adding both terms in Fig. 6. 

4.3 Stochastic counterterm: ~ A; 4 

Finally, we consider the addition of a stochastic counterterm (14) to the one-parameter two- 
loop formula from (4). Naively, we expect this term to contribute only at a much higher order 
in perturbation theory than the order at which we are working. However, there is subtlety in 
determining the expected size of this counterterm, as first pointed out in [17] in the context of 
biased tracers. An equivalent way of writing (14) is 

Pstoch(k) ~ (27t) 2 ( -• (24) 

\k M J n 


23 














- 2-loop EFT, c“ (1 ) = 0.53 (&nl/[2/j Mpc ]) 2 

- with stochastic term, c 2 (1 j = 0.54 (k^/Ylh Mpc -1 ]) 2 , c st0 ch = -35171 (&nl/[2/i Mpc -1 ]) 7 



Figure 10: The prediction of the EFTofLSS at two loops with the addition of the stochastic counterterm. 
We find that the fc rea ch is improved to about breach — 0.34 h Mpc -1 , with a theory error, in purple, that 
could decrease the reach to fc reac h ~ 0.21 /iMpc -1 . However, the magnitude and the sign of the required 
stochastic counterterm seem to be in conflict with the theoretical expectations, as explained in the main 
text. 


Here n is the number density of the objects that most contribute to the stochastic noise. The 
derivatives acting on this term are expected to be suppressed by the inverse length scale of the 
same objects, which we call &M- Now, for the dark matter power spectrum we expect both of these 
scales to be of order fc NL , which is why (14) is written with c stoc h expected to be an order one 
number. In particular, this is true for the part of the counterterm that is supposed to correct the 
perturbative loops, which have no information about the non-perturbative halos. However, given 
the high number of powers in which these scales appear, one should be careful, because the small 
difference in these scales from &nl can make quite a difference. For example, if very massive halos 
are the ones that contributes, 1/n is very large and ku is very small; vice-versa if it is small halos 
that contribute. Therefore, the value of c stoc h from (14) could indeed be very far from order one, 
and this would affect at which order in perturbation theory the term becomes relevant. 

The result of adding a stochastic counterterm to the prediction from (4) is shown in Fig. 10. 
We can see that adding the stochastic counterterm makes the EFT match the data up to k ~ 
0.28 h Mpc -1 , where the cosmic variance of the simulation is about 3 x 10 -3 . This is worse 
than the full UV-insensitive prediction, and not much better than even adding a single quadratic 
counterterm. The numerical values of the parameters from this fit are 

c 2 s{1) = 0.54 (fc NL /(2/rMp C - 1 )) 2 , c stoch =-3.5 x 10 4 (fc NL /(2/rMpc -1 )) 7 . (25) 

Although the magnitude of c st0 ch is comparable to what is required to fix the UV contribution 
to Po. loop, it is not compatible because it has the wrong sign. Indeed, the fitted value of c stoc h, 
negative in our case, is composed of a sum of UV and finite contributions. As seen in (18), the UV 
contribution requires a positive c s t oc h- Furthermore, as we will explain now, the finite contribution 
to c stoc h has to be positive as well, so that the overall sign of c st0 ch is expected to be positive. 
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Figure 11: Sketch of the nonlinear power spectrum in a toy model where the linear power spectrum has 
only short scale power. The prediction of the EFT as k —> 0 is 

(27r) 2 C sto ch(fc/fcNL) 4 fcNL ^ SinCe 

the power spectrum must be positive for all k, this implies that c at0 ch > 0 . 


Let us imagine a toy universe (shown schematically in Fig. 11) where the linear power spectrum 
is non-vanishing only over a small region between k\ < k < k 2 - At a time late enough so that 
/cnl < k\, the vanishing of Pu{k ) in the perturbative regime implies that the EFT prediction for 
the power spectrum is very simple: 

Ceft ~ (A) W + (A) TT-3 + • • • • (26) 

V fc NL/ K N L \ fc NL/ KnL 

Remarkably, in this toy universe the prediction at long wavelengths is entirely dominated by the 
stochastic contributions! Since the power spectrum is the expectation value of |(5*,| 2 , it must be 
positive for all k. By taking the limit k —> 0, we conclude that c st0 ch A 0. Notice that this argument 
requires that c stoc h is the finite contribution. The moment c st0 ch includes a renormalization of a 
loop, we cannot make this argument any longer. But, as we discussed, this is not the case at hand. 

We therefore conclude that we have no evidence of the necessity of adding a stochastic coun¬ 
terterm before the other counterterms in the UV-insensitive calculation 15 . We have tried to add 
the stochastic term after the inclusion of just the quadratic counterterm associated with ci, End¬ 
ing the same conclusions as when we add the stochastic term without the C\ term, ft is possible 
that the stochastic term should be added after both the c\ and c 4 terms are incuded, where the 
EFTofLSS stops fitting the nonlinear power spectrum due to lack of power. It is also possible that 
the stochastic counterterm might play a relevant role before the fc rea ch of the calculation with C\ 
and C4, so that, without the stochastic term, the additional terms would lead to a lack of power in 
the EFT prediction at a lower wavenumber. Though this is possible, we have no evidence of this 
from the fit to the power spectrum until the estimated /c reac h of the computation. But in either of 
these cases, it is likely that one should evaluate the three-loop contribution first. 

15 It is also reasonable to ask whether the coefficient of the stochastic term is so large that it should be included 
before any two-loop term, but after the one-loop terms from (2). When we fit such a formula to the data, however, we 
find that the fc-reach is not significantly improved over the one-parameter one-loop prediction, and furthermore, the 
sign of c s toch, which again is just given by the finite contribution, is negative, just as for the “two-loop+stochastic” 
prediction. 
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Figure 12: The prediction of the EFTofLSS for the matter power spectrum as a function of redshift. In 
black, we plot the EFTofLSS with three counterterms, with parameters c^, c\, and C 4 fit separately 
at each redshift. The darker grey band corresponds to an estimate of the theoretical error estimated 
by taking the value of which is 1 <t off from the best fit obtained at 0.75/cfi t . We also plot two- 
loop EFTofLSS predictions with different combinations of counterterms, and various other lower-order 
predictions. We see that the L’ re ach is higher and higher with the higher redshifts, and the gain with 
respect to SPT is very substantial at all redshifts. In is expected that the fc reac h as a function of redshift 
is a smooth function of z, once we take the theoretical error in account. 


5 Higher Redshifts 

We now proceed to study redshifts higher than z = 0. This will be useful to explore how much the 
breach of the EFTofLSS is improved as we move to higher redshifts, and also to explore the time 
dependence of the counterterms. 

5.1 Fits to the power spectrum 

The results of applying the same procedure that we described at z = 0 to higher redshift are 
given in Fig. 12. Figures of the values of the counterterms as a function of fc max , from which we 
determine t and the theoretical error, are provided in App. B. When we consider the calculation 
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done with the three relevant counterterms (with parameters c^, ci, and C 4 ), we clearly see that as 
we move to higher redshifts, the fc-reach is relevantly improved, to k ~ 0 . 6 hMpc -1 at z = 1 and 
k ~ 1.1 /rMpc -1 at z — 2. Based on our estimates of the theoretical error in the fits, the fc-reach 
of the prediction could potentially be as low as k ~ 0.26 hMpW 1 at z = 0, k ~ 0.4/rMpc -1 at 
z — 1, and as low as fc ~ 0.9/iMpc ^ 1 at * = 2. Nevertheless, the EFT provides a substantial 
gain with respect to other analytical techniques, such as SPT at two loops (also shown in Fig. 12). 
Such a gain becomes very important once we consider that the number of available modes scales 
as breach 3 - For example, at redshift z = 1, the gain in number of modes with respect to two-loop 
SPT is about 200. At z = 2, this same number is closer to 400 lf> . 

5.2 Additional checks of fitting procedure 

Given that we are fitting the power spectrum, albeit with very small error bars, with three pa¬ 
rameters, there is a certain concern that we might be overfitting, despite the fact that we have 
designed the fitting procedure described in Sec. 2.2 to minimize this possibility. We try to limit 
the possibility of overfitting by performing the following additional checks. First, as mentioned, 
we present a theoretical error on the prediction. Second, we have verified that the functional forms 
of the various counterterms do not cancel each other relevantly 1 ' . Third, we checked that if we 
try to fit the numerical data by setting P 2 -I 00 P = 0 , we are unable to match the numerical data as 
successfully as when we include iVioop- 

Finally, we check that the numerical values of the parameters we obtain are consistent with the 
size that is induced by the uncontrolled UV. This can be estimated by repeating the procedure used 
to obtain (18), with the only difference that the counterterms are estimated by fitting the difference 
of P 2 -I 00 P computed with cutoff A = 00 and a ^-dependent A(z) that roughly approximates what 
the nonlinear scale is expected to be as a function of redshift. In more detail, we choose the cutoff 
to be A(z) = 2 k Teach (z) and fit the counterterms to P^“(/c, z) — P 2 A ioop reach( ^(&, z) over the k-range 
0.05 /iMpc -1 — 0.5/c reac h(z) including a 0.3% error bar on the computation of the integrals. We 
obtain the following values for the counterterms for z = { 0 , 1 , 2 } when including the counterterms 

16 In particular, one can notice that SPT fails to match the data at such low wavenumbers that cosmic variance 
plays a relevant role in determining its fc rea ch- A more accurate estimate of its fc re ach can be obtained by noticing 
when SPT significantly deviates from the EFT predictions. 

1( In [16] instead it was noticed that by pushing k ren to higher values, there was some cancellation among two-loop 
diagrams and the W_io 0p counterterm. This cancellation would have disappeared if some slightly different value 
of was chosen, implying a lower fc reac h of the theory. See Sec. 3.2 for further discussion. 
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associated to ci, c 4 and c stoc h 18 


4 UV) (0 = {0,1,2}) = {-1.1,-0.27,-0.064} (fc NL /(2 hMpc -1 )) * 2 , (30) 

cf V) (^ = {0,1,2}) = {-5.3,-0.46,-0.050} (k NL /(2 /iMpc -1 )) 4 , 
cffi(z = {0,1, 2}) = {6.2 x 10 4 , 2.4 x 10 3 ,130} (fc NL /(2 h Mpc" 1 )) ? , 

and 

cj uv) (^ = {0,1,2}) = {-1.6,-0.47,-0.092} (fc NL /(2/iMpc -1 )) 2 , (31) 

ci UV) (^ = {0,1, 2}) = {-7.0, -0.74, -0.068} (k NL /(2 h MpW 1 )) 4 , 

when fitting only the terms associated to c 4 and c 4 . Applying the same method to Pi-\oop(k, z) 
(fitting from /c min = 0.005 hMpc -1 ) we have 

c’§ V \z = {0,1, 2}) = {-2.2, -0.28, -0.11} (fc NL /( 2 h Mpc ” 1 )) 2 . (32) 

These values are in general not very different from the coefficients that we find fitting the power 
spectrum numerical data in Table 2 and the values of c 4 and c 4 are quite independent of the presence 
of the stochastic term. Note that the values of the parameters are also quite independent from the 
choice of A, not varying more that a factor of 2-3 when considering the range A = 1 — 2 /c rea ch(^)- 
They are also stable under the change of the fitting range as one varies /c max = 0.2 — 0.9 fc reac h- One 
also notices that while at z = 1 the quadratic counterterm seems to bring most of the fc-gain, this 
is not the case at z = 0 and z — 2. Thus, along with the arguments about UV-sensitivity, the data 
themselves seem to indicate that the inclusion of all three counterterms is the most appropriate 
choice. 

In summary, we find these checks to be quite successful. We conclude that we find no strong 
indications that we are overfitting, even though we acknowledge that a better determination of 
the value of the counterterms by analyzing higher statistics or observables, as done for example 
in [9, 10, 20] 19 , or with direct measurement from small 77-body simulations [2], or by including 
higher order terms, would be helpful. 

18 If instead we choose A = 1 fc r each(-z) we obtain 

4 UV) (z = {0,1,2}) = {-2.8, -0.73, -0.18} (fc NL /(2 /tMpc" 1 )) 2 , (27) 

cf V) (z = {0,1,2}) = {-13, -1.2, -0.14} (fc NL /(2 /iMpc -1 )) 4 , 
c stodl( 2: = {0, b 2}) = {1.4 x 10 5 , 5.6 x 10 3 ,340} (/cnl/(2 h Mpc -1 ))" , 

and 

4 uv) (^ = {0,l,2}) = {-3.9,-l.l,-0.20}(fc NL /(2/iMpc- 1 )) 2 , (28) 

cf v) (z = {0,1,2}) = {-17,-1.7,-0.15} (fc NL /(2^Mpc- 1 )) 4 , 

when fitting only the terms associated to c\ and c 4 . For Pi_i 00 p(7’, z) we have 

c 5q V) (~ = {°> b 2}) = {-1.2, -0.11, -0.04} (fc NL /(2/ l M pC - 1 )) 2 . (29) 
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z = 0 

z = 1 

z = 2 


c 2 

U(D 

Cl 

c 4 

c 2 

c dU 

Cl 

c 4 

c 2 

C s(l) 

Cl 

c 4 

oniy c 2 (1) 

0.53 

X 

X 

0.20 

X 

X 

0.073 

X 

X 

only c 2 (1) & Ci 

0.48 

0.47 

X 

0.18 

0.23 

X 

0.066 

0.063 

X 

only c 2 (1) & c 4 

0.48 

X 

-2.55 

0.19 

X 

-0.29 

0.069 

X 

-0.033 

c 2 (1) & Ci & c 4 

0.48 

-0.74 

-6.41 

0.18 

0.22 

-0.014 

0.060 

0.15 

0.040 


Tabic 2: Table of the numerical values of the counterterms as a function of redshift z and for the various 
combinations that are studied in the paper. The units are (&nl/(2 h Mpc -1 )) 2 for and ci and 
(fc N L/(2/ i Mpc- 1 )) 4 for c 4 . 


5.3 Time-dependence of counterterms 

Since time translations are spontaneously broken in the universe, the time-dependence of the EFT 
parameters is unknown. The timescale of these coefficients is expected to be of order Hubble, so 
there should exist approximate functional forms related to this timescale. For each counterterm, 
we present a quasi-two-parameter fitting function that works reasonably well. 

We parametrize the time-dependence of each counterterm as the sum of two power laws in the 
following way 


c 2 s(1) (z) = A s D L (z) a ‘ +B s D l (z , (33) 

Ci(z) = Ai Di(z) ai + Bi Di(z)^ 1 , 

c 4 (z) = A 4 D 1 (z)° 14 + B 4 D 1 (z)^ 4 . 

The first power law, characterized by Ai and cq, represents the expected time-dependence of the 
counterterm induced from the cancellation of UV part of the loops, while the second power law, 
characterized by B\ and /%, is expected to be associated to the finite terms. For each of the three 
counterterms c 2 ,-^, ci, and c 4 , we fit these four coefficients. However, we call this a quasi-two 
parameter fit because we constrain the values of Ai and cq to lie close to the values obtained from 
fitting the time-dependence of the UV coefficients from (27) by a power law. The time-dependence 
of the parameters depends on which counterterms are included in the power spectrum fits, as this 
changes the reach of the theory, which in turn changes the cutoff A(z) used in the determination of 
c^/jx, Ci, and c 4 . Here we use the k -reach obtained by including both the counterterms associated 
to Ci and c 4 as we believe that they are the ones that need to be included. We fix the range that Ai 
can take to be between 0.33 and 3.0 times the best fit of the UV part of the parameters q, and ai 
to lie within 0.75 and 1.33 times the best fit. These values were chosen by analyzing the change in 


19 We notice that these additional ways to determine the value of the counterterms are possible because the 
counterterms of the EFTofLSS are terms that appear in some equations of motion and for which we know their 
origin in terms of UV degrees of freedom. This implies the fact that the same parameter appears in multiple 
observables or that one can measure them using directly dark matter particles. This is one of the characteristics 
because of which the EFTofLSS, being a theory and not a model, is more predictive than other approaches. 
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Figure 13: Time dependence of c\ and C 4 according to (33), where A* and a* are constrained to 

be close to the values obtained from fitting (27). 


the best fit values occurring when varying the ratios A.(z) / k vea , c h(z) from 1 to 3 and also from the 
change in their values when changing the fitting range from k = 0.05 hMpc -1 to 0.2 k Teac h(z) to 
k = 0.05 hMpc -1 to 0.9A(z). We believe that these represent reasonably well the uncertainty on 
the determination of the UV contribution to c^x, Ci, and C 4 . The results of the fits are presented 
in Fig. 13. 

Finally, we point out that the ^-dependent fitting functions from (33) should be understood 
as representing the typical scaling of the counterterms. If one were the use the numerical value 
of the parameters obtained from Figs. 13, one would find that the fc rea ch of the theory is smaller 
than when the counterterms are fitted at each redshift, because the fitting functions are not an 
exact match to the parameter values obtained from the separate fits. In Fig. 14, we show the 
results of the comparison to data when the fitting functions are used. It is expected that by using 
more general fitting functions, and combining with measurements of these parameters from several 
independent large-scale statistics, or from small-scale degrees of freedom in simulations (as in [ 2 ]), 
one can afford to use fewer parameters for the counterterms than if one were to have if the fit was 
performed independently at each redshift. Again, we leave a detailed study of this to future work. 


6 Conclusions 

In this paper, we have performed a high-precision comparison between the prediction for the dark 
matter power spectrum from EFTofLSS and nonlinear measurements from an Wbody simulation 
with a large box size and high number of particles, from the Dark Sky simulation set. The much 
higher precision of the numerical data allows us to better study the contribution of the counterterms 
than it was possible to do in previous studies. Starting at redshift z = 0, we have evaluated the 
two-loop prediction with one free parameter, as previously presented in e.g. [4], finding that it 
matches the data up to wavenumber k ~ 0.15 hMpW 1 to a precision of 0.3%. The k ieSLc h is smaller 
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^—2-loop EFT with c\ , and C 4 fit separately at each redshift 

^—2-loop EFT with c 1 , and C 4 given by fitting functions 



Figure 14: Comparison of the two-loop EFT power spectra when the three parameters are either fit 
separately at each redshift (solid black curves) or obtained from the fitting functions from (33) (dashed 
red curves). When the fitting functions are used, the fc reac h of the prediction is decreased slightly, but will 
still likely be acceptable for a range of applications. 


than what was previously presented in [4], where the error bar was taken to be ~2%, because of 
the much higher precision of the available numerical data and of a lower choice of k ien that reduces 
the theoretical error (and also removes the accidental cancellation between several two-loop terms) 
that was accidentally improving the fc-reach. 

We have tried to quantify how much the two-loop calculation with only one counterterm was 
sensitive to the contribution from the non-linear modes. We have done this by checking the 
contribution to the result from modes with wavenumber between 2 h Mpc -1 and infinity, which 
are not under perturbative control. We find that the contribution from these modes is non- 
negligible at the current level of precision, ft can be cancelled by turning on some additional 
counterterms that were not present in the first calculation of [4], effectively associated to one 
quadratic couuterterm and one four-derivative one. The fact that the counterterms have the 
power of making the calculation UV-insensitive represents a non-trivial consistency check of the 
EFTofLSS. The size of the prefactors ci and C 4 of these counterterms gives an indication of how 
much the short distance dynamics should contribute to the generation of these counterterms, with 
the expectation that the actual coefficient should not be much different than this. We therefore 
conclude that these terms should be included in a consistent two-loop calculation, resulting in the 
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k [h Mpc 1 ] 

Figure 15: The prediction of the EFTofLSS at linear, one-loop and two-loop levels. At one loop, only one 
counterterm, the so called speed of sound c^, is used, while at two loops two additional counterterms, 
respectively a sort of non-linear and an higher derivative speed of sound, ci and C4, are used. We notice 
the order by order improvement of the theory, and the remarkable fc reac h. 


following expression for the two-loop power spectrum: 


-pEFT-2-loop(&, z) — -PeFT- 1-loop (^5 z) + [D 1 (z)] 6 -P2-loop(^) — 2(27r)c^ 2 ) (%) - - ^ P\l{k) 


vNL 


]4 p(c a ) 

J - r l-loo 

l4 rj(quad, 1), 


c + 1 


Wc s (i)(^[Di{z)\ P^’ op (k) + (2vr) ^1 + 2 ^ + 2 ^ J [c a{1) (z)] [D^z)] -j—jPutk) 


k A 


+ {2n)c 1 (z)[D 1 (z)] A P^ \k) + 2(2n) 2 c 4 (z)[D 1 (z)] 2 — 4 P 11 (k) 

^NL 


(34) 


We find that the reach of this prediction at z = 0 is fc rea ch — 0.34 hMpc -1 , where the cosmic 
variance of the simulation is roughly one per mill (see Fig. 15). This is exquisite precision, and a 
remarkable success of the EFTofLSS. All former techniques fail at one per cent at k ~ 0.04 h Mpc -1 , 
which implies a huge number, potentially even three orders of magnitude, of additional modes that 
are under analytical control. We also consider other combinations of counterterms, in an effort to 
investigate if some of them give a negligible contribution. We find that adding a single higher- 
derivative or quadratic counterterm alone is not enough to reproduce the /c reac h of the full UV- 
insensitive prediction, and that including multiple quadratic counterterms can improve the fc reac h, 
but this is likely the result of overfitting. We discuss the inclusion of a stochastic term that would 
be more important than the terms we considered so far, and we argue that the overall magnitude 
of the coefficient that would be needed by the fit is too large compared to what is expected from 
the UV physics; we also find the sign to be inconsistent with theoretical considerations. 

We then perform the same study at higher redshift, with results shown in Fig. 12. We fold that 
the fc reac h of the theory grows at higher redshift, performing remarkably better than former analyt¬ 
ical techniques. We also fold evidence that inclusion of all of the three counterterms is necessary 
in order to maximize the k reac h of the theory at all redshifts. The exploration of the theory at all 
redshifts allows us to study the time-dependence of the counterterms. We fold that a reasonable fit 
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to the time-dependence of each of the three parameters can be obtained with a fonr-parameter fit, 
the sum of two power laws. One of these power laws can be highly constrained by imposing that 
its size and time-dependence is compatible to the size of the UV-induced counterterm, estimated 
by how much the calculation is sensitive to scale not under perturbative control, above 2fc reach (^). 
This constrains quite strongly two of the four parameters for the time dependent fit, so that we 
call this a quasi two-parameter fit. The resulting functional form of the counterterm is not precise 
enough to match the /c reac h found at each redshift when performing an independent fit. It however 
does give a good estimate of the size of the counterterms at each redshift. 

Given that we are determining three free parameters from the fit to the power spectrum at a 
single redshift, the risk overfitting is quite realistic. We have therefore performed several consis¬ 
tency checks on the calculations, as described in Sec. 5.2 ( 20 ). We find no strong evidence that 
we are overfitting, or that the fc rea ch of the theory has not been reasonably estimated. Ultimately, 
measurements of higher 77 ,-point functions [10], or measurement of the counterterms directly from 
the UV degrees of freedom [2], as well as the addition of the next order terms, in particular of the 
three-loop power spectrum, will also help in addressing the uncertainty associated to the possibility 
of artificially inflating the match of the theory to the data. 

Comparing the EFTofLSS to numerical data at this level of precision raises concerns about 
the accuracy of both theoretical calculations and results from simulations. On the side of the 
EFTofLSS, it is relatively simple to check the convergence of the calculations, as we have direct 
control over to them. One systematic mistake that is performed on the EFTofLSS is the ap¬ 
proximation of the time-integral in the perturbative expressions with corresponding factors of the 
growth factor D±. This is an approximate result when dark energy is present, and so it becomes 
better and better with increasing redshifts. Several studies (e.g. [39]) have verified that this ap¬ 
proximation is accurate to the 10 -3 level for the one-loop correction to the power spectrum at 
z — 0. It is not particularly hard to perform the calculation in the EFTofLSS with the exact time 
dependence. This was done at one-loop in the EFTofLSS in [2], and an extension at two loops has 
been in the planning for some time [40]. 

On the simulation side, it is not clear if currently available simulations reach the required 
numerical accuracy, or if they have even been tested to the required level. In fact, a recent 
reference [41] shows 0.6% difference between different numerical codes for k < 1/rMpc -1 (even 
though the authors commit to quoting just less than 1%), and differences of 3% for k < 10 hMpc -1 . 
This reveals that per mill precision is a very far goal. For example, systematics associated to the 
growth factor could be important, as different terms have different powers of the growth factor. 
Furthermore, in the largest of the Dark Sky simulations, we have found that, if taken at face 
value, the power spectrum measured from the initial snapshot at k ~ 0.3/tMpcU 1 is different by 

20 In particular, we have tried to estimate the theoretical uncertainty due to lack of computation of the higher 
order terms, which is quite large. We have represented this as a shaded band connecting our best curve with the 
curve obtained if we fit the data up to 0.75 fcfit(-z), even though this should be meant as a very rough estimate 
of the error. We have also performed several additional sanity checks on the counterterms, such as checking that 
there are no unjustified cancellations between the several terms, checking that if we were to remove P 2 -loop from 
the calculation, we would not be able to fit as well the data at all Us as we do when including LLioop, and finally 
checking that the size of the counterterms that we obtain from the fit is compatible with what is expected to be 
induced from UV-physics. 
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Figure 16: The prediction of the EFTofLSS at linear, one-loop and two-loop levels, with the same 
parameters as in Fig. 15, but zoomed in at low wavenumbers. Regardless of the /c reac h ; we can see the 
remarkable order-by-order convergence of the theory at low fc’s, probably beyond the precision of the 
numerical simulations and surely obtained with much less computational cost. 


about 1% from the input linear power spectrum, by more than what expected by the fact that 
the initial conditions have been evolved with 2LPT 21 . Such a mismatch, if true and important 
at low redshift, would be a concern for the accuracy of the comparison. However, it seems to 
us that, in order to improve the accuracy to below one percent, not only substantial work on 
the numerical codes needs to be done, but also, as pointed out in [41], the numerical cost of the 
computation may significantly increase due to the number of employed time steps and particles. 
Here in this paper we assumed the absence of systematic errors in the simulations, which clearly 
is an oversimplification. However, until numerical simulation data are provided with an estimate 
of such errors, we believe this is the approach that runs the smallest risk of enhancing the reach 
of the EFT. Clearly, understanding the size of the systematic errors of the simulations is beyond 
the scope of the paper. 

To be more quantitative, in this paper we have matched the EFTofLSS prediction to numerical 
data at the level of 10~ 3 , even though the data could be affected by systematics at the level of 
ICC 2 . A systematic effect of order 1% or less could have potentially very large consequences for 
the EFTofLSS. For example, if we were to allow for percent deviations between theory and data 
at low wavenumbers, we could match the data to higher wavenumbers with fewer parameters, as 
in the first two-loop calculation of [4]. However, if we assume that the largest systematic effects 
of simulations happen by incorrectly describing the short-distance dynamics, but in a way that 
conserves matter and momentum, then the difference in the prediction in the EFTofLSS for the 
data obtained with the correct or the incorrect dynamics is completely re-absorbed in a difference 
in the counterterms, which need to be measured in any event. It is with this hope that we use the 
numerical data by accounting only for their cosmic variance. 

21 We did not investigate this mismatch in detail, but it seems likely that is related to how the window function 
associated with the measurement process is deconvolved from the measurements (an effect that is probably in fact 
degenerate with the tree-level EFT counterterm at leading order). 
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We conclude with the following. Much attention has been paid to the /e-reach of the EFTofLSS, 
rather than on its accuracy at low wavenumbers. This is because it is impossible to check how 
accurate is the theory below the precision and the accuracy of the numerical data. However, 
regardless of the true h-reach of the theory, it is very likely that its predictions can achieve extremely 
high precision at long scales (e.g. k < 0.1 hMpc -1 ), probably much higher than that of numerical 
simulations, and surely in a less computationally demanding way (Fig. 16). This will eliminate the 
need to measure such large-scale observables from large simulations, allowing more computational 
time to be spent on smaller and more accurate simulations of nonlinear physics. The amount of 
progress that occurred in the last couple of years, since the emergence of the EFTofLSS, is an 
indication that the study of Large Scale Structure is rapidly becoming a high-precision science. 
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Appendix 


A Comments on quadratic terms in the effective stress 
tensor 


In this appendix, we provide further discussion about the quadratic terms in the effective stress 
tensor that we consider. These terms were first listed in Eq. (8), and we repeat them below for 
convenience: 


( dr) Pl l D (1 — <5) x 1 9* — s'j , d l [d 2 4>] 2 , d l \dW k 4>djd k <f\ , cfcP§ <9j<9 2 0 j . (35) 


As explained in [4], if we compute correlation functions of the matter overdensity S, one can use 
the bare velocity held equations, in which the counterterms appear only in the form ( dr) Pl l = 
(1 + 5) _1 <9jT u and are solely associated to r u . This is the origin of the factor of (1 — 6) in front of 
the list of counterterms. 

For clarity, we give the derivation of the last counterterm in (35), which is the less obvious. 
We are going to show that this term appears in c^rA Starting from Ty,-, we can include 


T, 


V 
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We can add and subtract other terms that appear in to make this contribution simpler. We 
can write 


di 2 (pdjV 3 

'-nw) 


= d 2 (j) 


djV 1 


\5ijdtv 1 
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didj(f) - ^ Sij&qb 
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The terms d 2 (j) — (didjCp — |<5 i:) <9 2 0) j and d 2 cj) ^ {-n{a)J) ~ ^)) are third order, 


and contribute to the power spectrum as k 2 Pu(k). The terms in d 2 (j) y — hj and d 2 (p (— |5jj<9 2 0) 
are degenerate with the second counterterm in ( d 2 cj) ) 2 . The leaves us with d 2 (pdidj(j). When we 
consider its contribution to djT 13 , we have two terms: 


djr 12 


D dj d 2 (pdi dj 4> + d 2 


jd 2 (j) = djd 


2 ' k * i d j (p+ -di(d 2 4>) 2 


(38) 


The hrst term is the fourth counterterm, while the second is degenerate with second counterterm. 
We thank Matias Zaldarriaga for discussions about this point. 


B Parameters at higher redshifts 

In this short appendix, we provide Figures 17 and 18, which show the value of the fit parameters 
as a function of fcg t at z = 1 and z — 2. 
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